require(lubridate)
require(dhlvm)

issue.labels <- c("防衛力強化", "敵基地攻撃", "辺野古移設", 
                  "消費税率下げ", "累進課税強化", "法人税率上げ", 
                  "治安優先", "処理水放出", "夫婦別姓", 
                  "同性婚", "日米同盟強化", "格差より経済", 
                  "国内産業保護", "原発廃止", "憲法改正")

#### データの準備 ####
## 東京大学谷口研究室・朝日新聞2021年衆院選政治家調査データ
elite.data <- read.csv("2021UTASP20211126_excludedQ11.csv")
# 対象者数（本文で言及）
nrow(elite.data)
# 回収率（本文で言及）
round(mean(elite.data$RESPONSE), 3)
# 全項目無回答者を除外
elite.data <- subset(elite.data, RESPONSE == 1)
# 回答者数（本文で言及）
nrow(elite.data)

## 本研究で収集した有権者調査データ
voter.data <- read.csv("main_data.csv")
# 回答者数（本文で言及）
nrow(voter.data)

# 学歴の分布（付録B.1）
voter.data$edu.cat <- ifelse(voter.data$edu < 3, 1,
                             ifelse(voter.data$edu < 6, 2, 3))
# 1: 低学歴，2: 中学歴，3: 高学歴
round(prop.table(table(voter.data$edu.cat)), 3)

# 比例投票先の分布（表A1）
# 人数
addmargins(table(voter.data$PR))
# 全体での割合
round(prop.table(table(voter.data$PR)), 2)
# DKを除いた割合
round(prop.table(table(voter.data$PR[voter.data$PR != 12])), 2)
# 投票者における割合
round(prop.table(table(voter.data$PR[voter.data$PR < 11])), 2)

# ロシアのウクライナ侵攻後の回答の割合（付録B.4）
voter.data$day <- day(voter.data$EndDate)
voter.data$hour <- hour(voter.data$EndDate)
voter.data$end.time <- voter.data$day * 100 + voter.data$hour
round(mean(voter.data$end.time < 2412), 3)

#### LDA-Sのパラメータの推定 ####
## LDA-Sに使うデータを抽出
# 候補者の争点態度
y1 <- cbind(elite.data$Q6_1, elite.data$Q6_2, elite.data$Q6_5, 
            elite.data$Q6_9, elite.data$Q6_10, elite.data$Q6_11, 
            elite.data$Q6_13, elite.data$Q6_14, elite.data$Q6_15, 
            elite.data$Q6_16, elite.data$Q7_1, elite.data$Q7_3, 
            elite.data$Q7_4, elite.data$Q7_5, elite.data$Q8)
# 欠損値を6に置換
y1[y1 == 99] <- 6

# 有権者調査の回答者の争点態度
y2 <- cbind(voter.data$P1_1, voter.data$P1_2, voter.data$P1_3,
            voter.data$P1_4, voter.data$P1_5, voter.data$P1_6,
            voter.data$P1_7, voter.data$P1_8, voter.data$P1_9,
            voter.data$P1_10, voter.data$P2_1, voter.data$P2_2,
            voter.data$P2_3, voter.data$P2_4, voter.data$const1)

# 有権者調査の回答者の争点態度の分布（付録B.4，図A2）
preference.prop <- apply(y2, 2, function(x) prop.table(table(x)))

cairo_pdf("Figure_A2.pdf", width = 5, height = 2.5, pointsize = 7, family = "Meiryo")
layout(matrix(1:15, 3, 5, byrow = TRUE))
par(mar = c(1, 2, 2, 1), lwd = 0.5)
for (i in 1:15) {
  plot(NULL, NULL, bty = "n", xlim = c(0.5, 6.5), ylim = c(0, 0.4), 
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = seq(0, 0.4, 0.1), lty = 3, col = "gray80")
  for (j in 1:6) {
    polygon(c(j - 0.5, j + 0.5, j + 0.5, j - 0.5), 
            c(0, 0, preference.prop[j, i], preference.prop[j, i]), 
            col = ifelse(j < 6, "gray", "white"))
  }
  title(issue.labels[i])
  axis(2, lwd = 0.5)
}
dev.off()

# 候補者と有権者調査の回答者のデータをつなげる
y <- rbind(y1, y2)

# 候補者の公認政党（有権者調査の回答者は別カテゴリーの12を割り当て）
group <- c(elite.data$PARTY, rep(12, nrow(voter.data)))
G <- max(group)

## K = 7, ..., 12として，それぞれ5通りの異なる初期値を設定してパラメータを推定
for (i in 7:12) {
  for (j in 1:5) {
    combined.result <- list()
    set.seed(100 * j + i)
    K <- i
    eta <- list()
    for (k in 1:ncol(y)) {
      eta[[k]] <- matrix(1, K, 6)
    }
    alpha <- matrix(1, G, K)
    combined.result <- hlcModel(y, group, eta, alpha, 10000, 2000, 1)
    save(combined.result, file = paste0("combined_result_", i, "_", j, ".Rdata"))
  }
}

# BICの比較
combined.BIC <- matrix(NA, 6, 5)
for (i in 7:12) {
  for (j in 1:5) {
    load(paste0("combined_result_", i, "_", j, ".Rdata"))
    posterior <- posteriorMeans(combined.result)
    combined.BIC[i - 6, j] <- bic(as.matrix(y), group, 
                                  posterior$pi, posterior$beta)
  }
}
round(combined.BIC)
which(combined.BIC == min(combined.BIC), arr.ind = TRUE)
# 3行3列の結果（K = 9の3番目の推定結果）を採用

#### 各政策クラスターの政策選好 ####
load("combined_result_9_3.Rdata")
# パラメータの事後平均
combined.posterior <- posteriorMeans(combined.result)
names(combined.posterior$beta) <- issue.labels

# 元の結果を見て，解釈しやすいように政策クラスターを並べ替える
table(group, combined.posterior$z_assign)
combined.posterior$z_assign.2 <- c(4, 5, 2, 8, 7, 1, 9, 6, 3)[combined.posterior$z_assign]

## 各政策クラスターのβの推定結果（表2，図A3）
cluster.order <- c(6, 3, 9, 1, 2, 8, 5, 4, 7)
cluster.label <- c("自民C", "立憲C", "公明・国民C", 
                   "左派政党C", "維新C", "右派有権者C", 
                   "中間C", "左派有権者C", "DKC")

cairo_pdf("Figure_A3.pdf", width = 5.5, height = 7, pointsize = 5, family = "Meiryo")
layout(matrix(1:(15 * 9), 15, 9, byrow = TRUE))
par(mar = c(0.2, 0.2, 0.2, 0.2), oma = c(0, 3, 3, 0), lwd = 0.5)
for (i in 1:15) {
  for (j in 1:9) {
    plot(NULL, NULL, bty = "n", xlim = c(0.5, 6.5), ylim = c(0, 1), 
         xlab = "", ylab = "", xaxt = "n", yaxt = "n")
    for (k in 1:6) {
      polygon(c(k - 0.5, k + 0.5, k + 0.5, k - 0.5), 
              c(0, 0, combined.posterior$beta[[i]][cluster.order[j], k], 
                combined.posterior$beta[[i]][cluster.order[j], k]), 
              col = ifelse(k < 6, "gray", "white"))
    }
  }
  mtext(issue.labels[i], side = 2, line = 0.9, 
        outer = TRUE, at = 1 + 1 / 30 - i / 15, cex = 1)
}
for (i in 1:9) {
  mtext(cluster.label[i], side = 3, line = 1, 
        outer = TRUE, at = i / 9 - 1 / 18, cex = 1, font = 2)
}
dev.off()

#### 政策クラスターの分類結果 ####
## 候補者の政策クラスターの分類結果
# 政治家調査データに政策クラスター所属を挿入する
elite.data$cluster <- factor(combined.posterior$z_assign.2[1:nrow(elite.data)], 
                             levels = 1:9)

# 候補者の公認政党と政策クラスターのクロス集計（表2，表A3）
addmargins(table(elite.data$PARTY, elite.data$cluster))
round(prop.table(addmargins(table(elite.data$PARTY, elite.data$cluster), 
                            margin = 1), margin = 1), 2)

## 有権者調査の回答者の政策クラスターの分類結果
# 有権者調査データに政策クラスター所属を挿入する
voter.data$cluster <- factor(combined.posterior$z_assign.2[(nrow(elite.data) + 1):
                                                             (nrow(elite.data) + nrow(voter.data))], 
                             levels = 1:9)

# 「覚えていない」人を除外した比例投票先の変数を作る
voter.data$PR.2 <- voter.data$PR
# 除外される人数
sum(voter.data$PR == 12)
voter.data$PR.2[voter.data$PR.2 == 12] <- NA

# 有権者調査の回答者の比例投票先と政策クラスターのクロス集計（表2，表A4）
addmargins(table(voter.data$PR.2, voter.data$cluster))
round(prop.table(addmargins(table(voter.data$PR.2, voter.data$cluster), 
                            margin = 1), margin = 1), 2)
# 政策クラスター側からみた割合（本文で言及）
round(prop.table(addmargins(table(voter.data$PR.2, voter.data$cluster), 
                            margin = 2), margin = 2), 2)

## クロス集計のモザイクプロット（図1）
# 小政党をその他にまとめる
voter.data$PR.plot <- ifelse(voter.data$PR.2 > 8 & voter.data$PR.2 < 11, 8, 
                             ifelse(voter.data$PR.2 == 11, 9, voter.data$PR.2))

# モザイクプロット
mosaicplot.data <- table(voter.data$PR.plot, voter.data$cluster)
rownames(mosaicplot.data) <- c("自", "立", "公", "共", "維", 
                               "国", "れ", "他", "棄")
colnames(mosaicplot.data) <- c("自民C", "立憲C", "公明・国民C", 
                               "左派政党C", "維新C", "右派有権者C", 
                               "中間C", "左派有権者C", "DKC")

cairo_pdf("Figure_1.pdf", 
          width = 4.2, height = 2.6, pointsize = 7, family = "Meiryo")
par(mar = c(0.5, 0.5, 0.5, 0), oma = c(0, 0.5, 1, 0), lwd = 0.5)
mosaicplot(mosaicplot.data, main = "", 
           color = c("gray50", "gray80", "white"), las = 1)
mtext("比例投票政党", side = 3, line = 0, outer = TRUE, cex = 0.8, font = 2)
mtext("政策クラスター", side = 2, line = -0.5, outer = TRUE, cex = 0.8, font = 2)
dev.off()

## 長期的党派性による分析結果（付録C.3.3）
# 有権者調査の回答者の長期的党派性と政策クラスターのクロス集計（表A6）
addmargins(table(voter.data$PID, voter.data$cluster))
round(prop.table(addmargins(table(voter.data$PID, voter.data$cluster), 
                            margin = 1), margin = 1), 2)

#### 政治家調査のみを用いた推定結果との比較 ####
G <- max(elite.data$PARTY)

## K = 2, ..., 10として，それぞれ3通りの異なる初期値を設定してパラメータを推定
for (i in 2:10) {
  for (j in 1:3) {
    elite.result <- list()
    set.seed(100 * j + i)
    K <- i
    eta <- list()
    for (k in 1:ncol(y1)) {
      eta[[k]] <- matrix(1, K, 6)
    }
    alpha <- matrix(1, G, K)
    elite.result <- hlcModel(y1, elite.data$PARTY, eta, alpha, 10000, 2000, 1)
    save(elite.result, file = paste0("elite_result_", i, "_", j, ".Rdata"))
  }
}

# BICの比較
elite.BIC <- matrix(NA, 9, 3)
for (i in 2:10) {
  for (j in 1:3) {
    load(paste0("elite_result_", i, "_", j, ".Rdata"))
    posterior <- posteriorMeans(elite.result)
    elite.BIC[i - 1, j] <- bic(as.matrix(y1), elite.data$PARTY, 
                               posterior$pi, posterior$beta)
  }
}
round(elite.BIC)
which(elite.BIC == min(elite.BIC), arr.ind = TRUE)
# 6行1列の結果（K = 7の1番目の推定結果）を採用

## 政治家調査のみを用いた政策クラスターの分類結果
load("elite_result_7_1.Rdata")
# パラメータの事後平均
elite.posterior <- posteriorMeans(elite.result)
names(elite.posterior$beta) <- issue.labels

# 元の結果を見て，解釈しやすいように政策クラスターを並べ替える
table(elite.data$PARTY, elite.posterior$z_assign)
elite.posterior$z_assign.2 <- c(7, 4, 3, 5, 2, 1, 6)[elite.posterior$z_assign]

# 2つの推定結果との比較（表A5）
addmargins(table(elite.posterior$z_assign.2, 
                 combined.posterior$z_assign.2[1:nrow(elite.data)]))